/*
 * Programacao concorrente
 * Grupo-03a
 * Brenno Candido					7696500 
 * Jose Victor Uliana Martins		7656620
 * Maria Fernanda Garbim Avelino 	7277562
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <time.h>

#include "rand_bm.c"

#define NTHREADS 10

pthread_t thread[NTHREADS];
pthread_mutex_t mutex;

double S, E, r, sigma, T; /** Varias de entrada **/
double mean, std_dev; /** Variaveis que serao concorridas pelas threads **/
double *trials;
int M; /** Variavel de entrada: tamanho do vetor trials **/

/*
 * Standard deviation (desvio padão)
 * Esta função é usada para calculo do desvio padrão,
 * porem ela não a calcula diretente, ela apenas calcula
 * a soma dos quadrados de trials[i]-mean.
 * A parte final na funcao e concorrida por varias threads,
 * por essa razão o uso do multex, somando na varials globals
 * que no main ira finalizar o calcudo do devio padrão total
 * de trials.
 */
void* stddev(void *arg) {
	double sum = 0;
	int i, init, end; /** init e end são os limites que cada thread pode acessar de trials para não haver conflitos **/
	int nthr = (int) arg;
	
	/** calculo dos limites do loop **/
	init = nthr*M/NTHREADS; // inicio
	end = init+(M/NTHREADS); // fim

	for(i=init; i<end; i++) {
		sum += pow(trials[i]-mean, 2);
	}
	
	/** zona concorrida pelas threads **/
	pthread_mutex_lock(&mutex);
	std_dev += sum;
	pthread_mutex_unlock(&mutex);
	
	pthread_exit(NULL);
}

/*
 * Funcao que realiza o algoritmo do black scholes utilizando
 * o metodo de monte carlo.
 * Esta funcao ira rodar em thread preenchendo o vetor global
 * trials alem de calcular de ajudar no calculo da madia do
 * vetor trials.
 * Cada thread sera responsavel por uma porcao do vetor trials
 * para que ele posso ser preenchido em paralelo.
 */
void* BlackScholes(void* arg) {
	struct BoxMullerState state; //responsavel por gerar numeros aleatorios do metodo
	double t, aux, sum = 0;
	int i, init, end; /** init e end serao responsaveis pelos limite do loop **/
	int nthr = (int)arg; // numere da thread
	
	initBoxMullerState(&state); // inicializa a geracao dos numeros randomicos

	/** calcula so limites de cada thread **/
	init = nthr*M/NTHREADS;
	end = init+(M/NTHREADS);
	for(i=init; i<end; i++) {
		t = S*exp( (r-0.5*sigma*sigma)*T + (sigma*sqrt(T)*boxMullerRandom(&state)));
		if(t-E > 0) {
			aux = (-1)*r*T;
			aux = exp(aux);
			trials[i] = aux*(t-E);
		}
		else {
			trials[i] = 0;
		}
		// soma de trials para ser tirado a media
		// cada thread tera sua propria soma
		sum += trials[i]; 
	}

	/** zona concorrida pelas threads **/
	pthread_mutex_lock(&mutex);
	mean += sum;
	pthread_mutex_unlock(&mutex);
	
	pthread_exit(NULL);
} 

int main() {
	double confwidth, confmin, confmax;
	int i;

	// leitura da entrada
	scanf("%lf %lf %lf %lf %lf %d", &S, &E, &r, &sigma, &T, &M);
	
	printf("S\t%.4lf\n", S);
	printf("E\t%.4lf\n", E);
	printf("r\t%.4lf\n", r);
	printf("sigma\t%.4lf\n", sigma);
	printf("T\t%.4lf\n", T);
	printf("M\t%d\n", M);
	
	// alocacao de trials
	trials = (double*)malloc(sizeof(double)*M);
	if(!trials) {
		printf( "ERROR\n");
		return -1;
	}


	mean = 0;
	pthread_mutex_init(&mutex, NULL); // inicializa multex
	for(i=0; i<NTHREADS; i++) {
		pthread_create(&thread[i], NULL, BlackScholes, (void*)i); // lanca as threads
	}

	for(i=0; i<NTHREADS; i++) {
		pthread_join(thread[i], NULL); // espera as threads acabarem de serem executadas
	}

	// calcula a media de trials
	mean = mean/M;

	std_dev = 0;
	for(i=0; i<NTHREADS; i++) {
		pthread_create(&thread[i], NULL, stddev, (void*)i); // lanca as threads
	}

	for(i=0; i<NTHREADS; i++) {
		pthread_join(thread[i], NULL); // espera as threads acabarem de serem executadas
	}
	
	// calcula desvio padrao de trials
	std_dev = sqrt(std_dev/M);

	/** calculo do intervalo de confianca **/
	confwidth = 1.96*std_dev/sqrt(M);
	confmin = mean - confwidth;
	confmax = mean + confwidth;

	printf("Confidence interval: (%.4lf, %.4lf)\n", confmin, confmax);

	return 0;
}
